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Abstract 

This paper is concerned with linear quantum networks of N nondegenerate opti¬ 
cal parametric amplifiers (NOPAs), with N up to 6, which are interconnected in a 
coherent feedback chain. Each network connects two communicating parties (Alice 
and Bob) over two transmission channels. In previous work we have shown that 
a dual-NOPA coherent feedback network generates better Einstein-Podolsky-Rosen 
(EPR) entanglement (i.e., more two-mode squeezing) between its two outgoing Gaus¬ 
sian fields than a single NOPA, when the same total pump power is consumed and 
the systems undergo the same transmission losses over the same distance. This paper 
aims to analyze stability, EPR entanglement between two outgoing fields of interest, 
and bipartite entanglement of two-mode Gaussian states of cavity modes of the N- 
NOPA networks under the effect of transmission and amplification losses, as well 
as time delays. It is numerically shown that, in the absence of losses and delays, 
the network with more NOPAs in the chain requires less total pump power to gen¬ 
erate the same degree of EPR entanglement. Moreover, we report on the internal 
entanglement synchronization that occurs in the steady state between certain pairs 
of Gaussian oscillator modes inside the NOPA cavities of the networks. 


1 Introduction 

Entanglement is a key factor in certain quantum information processing tasks, such as 
quantum teleportation DEI- Continuous-variable quantum information is highly moti¬ 
vated, as preparation of continuous-variable entangled states is efficient and mathematical 
description of a continuous variable system is adapted to various physical systems (e.g., 
quadrature operators of light and total angular momentum operators of an ensemble of 
atoms satisfy the canonical commutation relations) [31I11E]- In particular, Gaussian states 
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play an important role in continnons-variable qnantnm information ascribable to 1) phys¬ 
ically, the gronnd state of a qnantized electromagnetic held is a Ganssian state; 2) mathe¬ 
matically, thongh an arbitrary Ganssian state lives in an inhnite Hilbert space, it is simply 
and completely characterized by the mean and (symmetrized) covariance of held operators, 
which are represented by a vector and a matrix with hnite dimensions, respectively; 3) ex¬ 
perimentally, entanglement in Ganssian states is easily obtained [SIEIEIIH]. For instance, 
Ganssian entangled beams can be generated by employing a strong coherent pnmp beam 
to shine a nonlinear crystal inside the cavity of a nondegenerate optical parametric 
ampliher (NOPA). Interactions between the pnmp beam and two modes in the cavity pro- 
dnce a pair of entangled ontgoing helds in Ganssian states which have strong correlations 
in qnadratnre-phase amplitndes [8]. Snch entanglement is called Einstein-Podolsky-Rosen 
(EPR) entanglement. 

In this paper, we investigate EPR entanglement between two propagating continnons- 
mode light helds and bipartite entanglement of two-mode Ganssian states. The former one, 
namely, the two-mode sqneezing of the two continnons-mode helds, is characterized by a 
two-mode sqneezing spectrnm V{iu:), which can be obtained via an associated qnantnm 
Langevin eqnation. Strong EPR entanglement is indicated by a high degree of two-mode 
sqneezing over a wide freqnency range. The snfhcient condition of EPR entanglement 
existing at a certain freqnency for a pair of qnadratnre-phase amplitndes is that the valne 
of the corresponding V{tu) satishes a snm criterion [HI E] . On the other hand, entanglement 
of two-mode Ganssian states can be measured by the logarithmic negativity E{t) in the 
time domain. The logarithmic negativity is calculated based on the covariance matrix 
of the position and momentum operators corresponding to the two modes. For strong 
entanglement, E(t) has a high positive value [Zliin]. More details of the two quantities of 
entanglement are given in Section 

In reality, a quantum system is dissipative due to interactions between the system and 
its environment. Such inevitable losses degrade entanglement. Thus, transmission distance 
and communication quality are limited. Failure of the communication may even happen 
nmn. Gonsequently, reliable generation and distribution of entanglement in the presence 
of losses in transmission channels is a central issue in quantum communications. 

Our previous work [12] has investigated EPR entanglement prepared by a dual-NOPA 
coherent feedback network where the two NOPAs are distributed separately at two distant 
communicating ends (Alice and Bob), as shown in Fig. (b). The degree of EPR entan¬ 
glement is influenced by the total pump power applied to the system, values of the damp¬ 
ing rates of the NOPAs, amplihcation losses induced by unwanted interactions between a 
NOPA and its environment, as well as transmission losses caused by leakage of photons 
along the transmission channels. Moreover, time delays in the process of transmission nar¬ 
row the bandwidth of suppressed two-mode squeezing. Gompared to a single NOPA placed 
in between the two ends (at Gharlie’s) over the same transmission distance as indicated in 
Fig. a (a), the coherent feedback conhguration consumes less total pump power to achieve 
the same degree of EPR entanglement when amplihcation losses are ignored and the damp¬ 
ing rates of the systems are identical; on the other hand, the dual-NOPA system proposed 
in [121 achieves higher level of EPR entanglement against transmission losses under the 
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Figure 1: (a) Entangled pairs (blue circles) generated by a single NOPA which is placed 
in between Alice (A) and Bob (B), say at Charlie’s (C). (b) Entangled pairs generated 
by a dual-NOPA coherent feedback system, in which two NOPAs are deployed at Alice’s 
and Bob’s separately, (c) Entangled pairs generated by an A^-NOPA coherent feedback 
system, in which N NOPAs are evenly distributed in a linear fashion between Alice and 
Bob. Transmission losses of the three systems are denoted by beamsplitters (BS) with 
transmission rates and a, respectively. (Calculations of transmission rates 

are given in Section]^) 
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same total pump power and damping rates. More precisely, for sufficiently high transmis¬ 
sion rates one would employ a distributed version of the scheme, while for low transmission 
rates a centralised architecture would be employed to give enhanced entanglement under 
the same pump power and damping rates. 

This paper studies a coherent feedback conhguration of up to six NOPAs that are evenly 
deployed between two ends (Alice and Bob) and connected in a linear coherent feedback 
interconnection, as shown in Fig. (c). Descriptions of optical components employed in 
the system as well as the dynamics of the system under the effects of losses and time delays 
are given in Section Section gives an analysis of stability conditions in systems without 
and with losses. We prove a theorem showing that the stability thresholds can be obtained 
by solving certain polynomials. Sectionj^is devoted to investigating and comparing degrees 
of end-to-end EPR entanglement of the A^-NOPA systems in absence of time delays but 
with losses considered. We give the ideal values of 9a and 9b as shown in Fig for WNOPA 
systems with the even number and the odd number of NOPAs. Moreover, we look into 
entanglement of two-mode Gaussian states related to pairs of optical cavity modes. In 
this section, there are some qualitative Endings such as in the absence of losses, the degree 
of end-to-end entanglement increases as the amplitude of pump beam for each NOPA 
approaches the value at which the system just loses stability; in the presence of losses, there 
exists an optimal value of pump amplitude of NOPAs for the end-to-end entanglement; 
the entanglement generated between collective modes is difierent between systems with an 
even and odd number of NOPAs; with the same pump power and without losses, increasing 
the number of NOPAs improves the end-to-end entanglement but not the entanglement 
of internal two-mode Gaussian states. The interesting qualitative observations described 
above are not straightforward to analyze quantitatively, so their quantitative analyses are 
left as topics for future research. Section explores efiects of time delays on the stability 
and entanglement of the systems. Finally, we restate the main results of this paper as a 
conclusion in Section |6l 

The following notations are adopted in this paper: i denotes \/^, the transpose of 
a matrix of numbers or operators is denoted by (-j^, and (■)* denotes (i) the complex 
conjugate of a number, (ii) the conjugate transpose of a matrix, as well as (hi) the adjoint 
of an operator. denotes an n by n identity matrix, Omxn is an m by n zero matrix (we 
simply write Om, if m = n), trace operator is represented by Tr[-], 5ij denotes the Kronecker 
delta and 5{t) denotes the Dirac delta function. (■) denotes quantum expectation and a(-) 
denotes the largest singular value of a matrix. 


2 System Model 

In this section, a brief introduction of quantum optical devices and dynamics of the system 
is given. Fig. shows a detailed system model of our A^-NOPA coherent feedback network. 
We take account of the system undergoing losses and time delays. The first NOPA (Gi) and 
the A^-th NOPA (Gjv) are placed at Alice’s and Bob’s, respectively; the other NOPAs are 
deployed in a linear line in between Alice and Bob. The length of the path between every 
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Figure 2: An A^-NOPA coherent feedback system undergoing losses and time delays. 


two neighbouring NOPAs is identical. Transmission losses are modelled by beamsplitters. 
The time delay in each path is denoted by a constant r. Two adjustable phase shifters 
with phase shifts 9a and 9b, respectively, are placed at two outputs separately in order to 
obtain the best two-mode squeezing between helds ^out,a and ^out,b- 


2.1 Quantum optical components 

2.1.1 NOPA 


Fig. I (a) gives a block diagram representation of a NOPA (Gj). Four ingoing boson 
helds are in the vacuum state, among which ^ioss,a,i and ^ioss,b,i are unwanted amplihcation 
losses. The held operators comply with the commutation relations, that is, for a boson 
held operator we have [ 6 (f)) 0('®)*] ~ ~ s). The main component of NOPA is a 

two-ended cavity with two orthogonally polarized boson modes a* and 6 , which obey the 
commutation relations [a,, a*] = 6 ij, [bi,b*] = 6 ij, [ai,b*] = 0 and [ai,bj] = 0. The helds 
6 n,a,i, 6 n,fe,i, ^ioss,a,i and interact with the modes via coupling operators Li = ^/yai, 

L 2 = y/lh, L 3 = \/~nai and L 4 = ^/i^bi, respectively, where 7 and k are damping rates. To 
yield EPR entanglement between outgoing Gaussian helds ^out,a,i and iout,b,i, a strong pump 
beam in a coherent state is employed to shine the nonlinear crystal inside the cavity. 
The modes are coupled with the beam via the system Hamiltonian H^ys = {^Ib* — aibi), 

where e is a real parameter related to the ehective amplitude of the pumping held [HI fT3] . 
The time-varying interaction Hamiltonian between the system and its environment is 

hfint(f) ^( 6 )T.(f) ^ ^ where \^in,a,i(t), ^in,b,i(t'), ^loss,a,i(t'), ^loss,b,i(t')^ and 

L = [Li, L 2 , L 3 , L^]'^ [13]. The time evolution of the mode and held operators in the 
Heisenberg picture are given by [T3| 


ttiit) = U{t)*aiU{t), 

bi{t) = u{tyhu{t), 
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where U{t) = ex^ {—i Hi^t{s)ds) is a unitary process satisfying the quantum white 
noise equation U{t) = {t) . Using the rules of quantum stochastic calculus, the 

dynamics of the NOPA Gi is described by the quantum Langevin equations PUSlUSlES] 

di (t) = - (^) ai (t) + (t) - {t) - \fKiioss,a,i {t) , 

h{t)= - , ( 1 ) 

and the output helds are 

^out,a,i {t') (^) A ^in,a,i (^) ) 

Lut,b,i (t) = (^) + iin,b,i {t) ■ ( 2 ) 


^in,a,i ^ 


^OUtrCri f 





^loss,a,i ^ 

^loss,b,i ^ 



r 

> 

\ 

f 


(a) (b) 

Figure 3: (a) A NOPA. (b) A beamsplitter. 


2.1.2 Beamsplitter 

The transmission loss in a path between two adjacent NOPAs is caused by loss of photons. 
It is modelled by a beamsplitter with transmission rate a and reflection rate /? = \/l — 
[isiin]. As shown in Fig. |^(b), output signal ^BS,out of the beamsplitter is the combination 
of the two ingoing fields, that is, ^BS,out = In our case, ^Bs,i is a white-noise 

held operator in the ground state and ^BS,in is an outgoing Gaussian held of a NOPA. 

Assume the total transmission distance between Alice and Bob is d kilometres. There¬ 
fore, the length of each path between every two adjacent NOPAs is km. Based on the 
assumption of that there exists around 0.2 dB per kilometre transmission loss in optical 

_ O.Old 

hbre [18], we obtain that the transmission rate of each beamsplitter is a = 10 . 

2.2 The N-NOPA coherent feedback network 

Based on the dynamics of a NOPA and the transformation of a beamsplitter as given above, 
an A^-NOPA coherent feedback system undergoing transmission losses and time delays as 
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depicted in Fig. |^has the following dynamics at time t > {N — l)r, when all the nodes of 
the network are connected, 


a'l (t) 
di (t) 


bN (t) 


bj (t) 


with outputs 


^out,b (t) 


Cout,a (t) 


^ 2 ^ (^) \/7'Cm,a,l (t) '\/~i^^loss,a,l (t) j 


(^) + (^) “ V^iloss,a,i (t) 


i— 1 

-7 ^ a’'ai-k{t - kr) - - {i - l)r) 

k=l 

2—1 

-ldy/l'^(y^~^iBS,a,i-k{t - kr), 
k=l 

^ ^ (t) + 2^^ '\/^Cin,h,N (t) 'y~^^loss,b,N (t) ; 

^ ^ bj (t) + 2®j (^) ~ V~^Closs,b,j (t) 

N-j 

-7 ^ - fer) - V7?m,fe,w(t - {N - j)T) 


k=l 


N-j 


^^BSAj+kit-kr), 


k =0 


N 


I - 1)'^) + «^“ ^mAiV {t- {N - l)r) 

V k=l 

N-l \ 

+/3 ^2 I , 

fc=l / 

/ iV 

e*^“ ( ^ak{t - {N - k)T) + {t - {N - l)r 

V fe=i 

Af-l \ 

+/3 ^ a^~^iBS,a,N-k{t -kr)] , 

k=l J 


(3) 


(4) 


where 1 < i < and 1 < j < iV. 

As a linear stochastic model, the dynamics of the system can be described by a linear 
quantum stochastic differential equation in the quadrature operators of the system [iniEo]. 
Note that quadratures of a bosonic mode, say Oj, are = a* + a* and af = — zoj + la*. 
Similarly, the quadratures of a held operator, say are = 6 + C ^^*4 = —i^i + iQ. 
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Define the vectors of quadratures corresponding to the iV-NOPA system as 

_ \ q nP nP If ... n'? nP }f F 

^ ~ L“l) '^1’^^2’^^2’"2 j'^2 5 ) “TV"’’ 

f = IfJ fP F’i pP pi pP 

^ ’ ^in^b^N ’ ’ ^/oss,a,l ’ ^/oss,a,l ’ 

tP d9 tP c<? tP 

^loss,b,l^ ‘»Zoss,6,l’ ^Zoss,a,2’ ^loss,a,2^ ^Zoss,b,2’ ‘’/oss,6,2’ ’ 

cP t9 tP tP 

^loss,a,N ’ ^loss,a,N'’ ^loss,b,N'’ ^loss,b,N'’ ^BS,a,l ’ ^i3S,a,l ’ 

t? tP tP ... t9 cP 

^BS,a,'2^ SB5,a,2’ ^BS,b,2^ ^BS,b,2^ > SB5,a,V-l> SB5,a,Af-l’ 

pq tP tq irP ]T 

^BS,b,N-l^ ^BS,b,N-l^ ^BS,b,N^ ^BS,b,N\ ’ 

[?oiii,a’?oiiZ,a> ?oiiZ,6’' 

In the absence of time delays, the system dynamics is of the form 


i (f) = Anz (f) + (f), 

(t) = Cnz (t) + Dn^ (t). 


(5) 

( 6 ) 
(7) 


For the A^-NOPA system, the covariance matrix of its 2A^-mode Gaussian state is 

PN{t) = ^Tr (^p(O) (^z{t)z{tf + {z{t)z{tf)^^'^ , Pn{0) = Po, ( 8 ) 

where p(0) is the initial density operator of the system. Furthermore, when time delays 
are neglected, PN{t) satishes the Lyapunov matrix differential equation 

= A^P^it) + PN{t)Al + BmBI. (9) 

Specially, the steady-state covariance matrix P^r = hmt^ooPv(^) satishes the Lyapunov 

equation mm 

AnPn + PnA'^ P BjsfBjf = 0 . ( 10 ) 


Dehne the parameters of the system as follows. For each NOPA, we set e = x'jr Hz and 
7 = y Hz, where x (0 < a; < 1) and y {0 < y < 1) are real parameters, and = 7.2 x 10^ 
Hz is a reference value for the transmissivity mirrors of the NOPA. Following [121 HI], we 
assume that k, = when e = O.Oy^ and the value of k, is proportional to the absolute 

value of e, so we set k = mm- Suppose that the total transmission distance is 

0.01 

1 km, then the transmission rate of each beamsplitter is a = 10 and time delay r in 
the path between any two adjacent NOPAs is around The range of the adjustable 

phase shifts 9a and 9b is (—7r,7r]. 


3 Stability analysis 


This section is devoted to the stability analysis of our A^-NOPA system. If the system 
is unstable, the mean total photon number in the cavity modes is continuously growing. 






which is undesirable. By stability, we mean that the mean of the quadrature vector z{t) 


becomes a zero vector as time approaches infinity, namely, {z{oo)) = 0, and Eq. (10) has 
a unique solution [21]. The system is stable when Apf in Eq. (|^ is Hurwitz, that is, all 
the eigenvalues of Ajv have real negative parts. We are interested in the range of x over 
which stability is assured. However, unlike the dual-NOPA system studied in our previous 
work na, it is infeasible to obtain an explicit expression for the stability threshold Xth at 
which the system just loses stability by checking the Hurwitz property of Aj^. Here, by 
regarding x as an uncertainty, we employ the p-analysis method from control theory 
to obtain the following lemma. 


Lemma 1 In the absence of time delays, an N-NOPA coherent feedback system is stable 
if and only if det{zojI — Hjv(x)) ^ 0 Vu G TZ and Vx G [0, Xth), where 0 < Xth < 1- 

Proof. First, for convenience, let us dehne the following matrices, 

0 


An — 


7 


K 


hN, Ao(x) — Jr 


0 

0 


0 

0 

0 


0 -- 
^ 2 


2 

0 

0 

0 


0 

0 


Ah = 


■ 0 

0 

0 

0 


-7 

0 

0 

0 ■ 

0 

0 

0 

0 

= 

0 

-7 

0 

0 

0 

0 

-7 

0 

0 

0 

0 

0 

0 

0 

0 

_ 


0 

0 

0 

0 


( 11 ) 


Recall that e = xjr, then from Eq. ([^ and Eq. ([^ the matrix An{x) of an iV-NOPA 
system is a 4A^ x 4A^ real matrix given by 


An{x) = 


Aq 

aAh 

a^Ah 


a^-\ 

aAa 

Ao 

aAb 


a^-f 

.N-2a 

■^a 


aAa 

Ao 

aA\ 

.N-l A 
■^a 


o?Aa 

aAa 

^0 


Ab 


+ AAr(x), 


( 12 ) 


where An{x) = /at ® Ao(x). To proceed, we now investigate the robust stability condition 
of a certain feedback system, as shown in FigQ The state space of the system is given by 

5 = A]^z + u, 

y = z, (13) 

where = A]^{x) — Aj^{x) and 5 is the state vector. Hence, the transfer function of the 
system is G{s) = {si — An)~^. As long as the stability condition of this system holds, our 
A-NOPA coherent feedback system is stable. 

According to Theorem 11.8 in [22], if the system (13) with transfer function G{s) is 


stable, then the closed-loop system in Fig. |^with structured uncertainty Ajv(x) satisfying 
||Ajv(x)||oo < 1 {v > 0) is internally stable if and only if 


sup ^{G{zu))) < r], 


(14) 
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■> 


Aj^ (x) 



u 




\/ 


Figure 4: A feedback system with transfer function G{s) and uncertainty AnIx). 


where 


li{G{iuj)) = 


min{o-(Ajv(3:))} 

0 , 


fordet(/ — G{iu:)Ai<[{x)) = 0, Vx € (0,Xt/j) 
fordet(/ — G{iu:)Am{x)) 7 ^ 0 , Vx G ( 0 ,xt/i). 


(15) 


First, we check the stability of the system (13). Using (11) and (12), we obtain that 
det(A/ —Am) = Hence G is stable because det(A/ — Am) has all its zeros in the 

left half plane by virtue of the fact that 7 > 0 and k > 0. Noting that 0 < x < 1, we have 

iN 

||Am(x)||oo = = I < 1^, therefore, V = Since ^ 

TZ, for (14) to be fulfilled we must have that det(/— G(*a;)AM(x)) 7 ^ 0 Vx G {0,Xth)- Since 
det(J — G(?C(;)Am(x)) = det(G(«a;)) det(za;J — Am(x)) and det(G(«a;)) = {ttu + ^ 

we obtain Lemma [H 


h^ivU)) 


= - > rjWu E 


Lemma shows that the system is stable as long as Am(x) does not have any purely 
imaginary eigenvalues. In fact, we shall prove that Am(x) has no eigenvalues on the imag¬ 
inary axis, which leads to the following theorem. 


Theorem 2 In the absence of time delays, an N-NOPA coherent feedback system is stable 
if and only if 0 < x < Xth, where Xth < 1 is the smallest positive root of the polynomial 
det (Am(x)). 


Proof. For convenience, we define m = n(x) 


^ and an invertible matrix 


Jm ® [ 1 0 0 0 ] 

Jm ® [ 0 0 1 0 ] 

Jm ® [ 0 1 0 0 ] 

Jm ® [ 0 0 0 1 ] 


(16) 


10 


















Exploiting (12), the characteristic polynomial of A]s[{x) is 


Pc{X, x) = det (A/4Ar - A 7 v(a;)) 


= det 


Al„(A) = 


A(A) = 


(L (A/4Ar — A]\f 

{x))L- 

‘) 



Al„(A) - 

n{x)lN 

On 

On 


-n{x)lN 

Alz(A) 

On 

On 


On 

On 

^n(A) 

n{x)lN 

On 

On 

n{x)lN Ai{X) 

N matrices given by 




A -I- m 0 

0 


0 


ay A -|- m 

0 


0 


a'^“^7 

ay 

X + m 

0 


a'^“^7 

a^y 

ay 

A -|- m 


A -|- m ay 

a^y 


a^“^y 


0 A -I- m 

ay 


a^“^y 


0 

0 

A -|- m 

ay 


0 

0 

0 

A -|- m 



( 17 ) 


(18) 


Since det 


A B 
C D 
DC and D is inver 


= det(AC' — BD) for any square matrices A, B, C, D such that CD = 
ible, we obtain 

2 


Pc(A, x) = (det (A„(A)y4z(A) - n{xflN)) , 
and 74„(A)74i(A) — n{x)‘^Iisf is a symmetric matrix given by 


(19) 


A„(A)A;(A) - n{x)'^lN 

(A + mp — n{xp a'y{\ + m) 

a7(A + m) {\ + mp - n{xp + l2,2 

q^ 7 (A + m) a (q! 7 (A + m) + (2,2) 


o^7(A + m) 
a (a7(A + m) + Z2,2) 
(A + mp — n{x)^ + Z3,3 


a^“^7(A + m) 
ia'y{\ + m) + 12,2) 
(a7(A + m) + [ 3 ^ 3 ) 


^7(A + m) ^ (a7(A + m) + ^ 2 , 2 ) (a7(A + m) + Z3,3) 

_ a-'^“'-7(A + m) a^~'^ laj{X + m) + 12,2) (a^lx + m) + [3^3) 


n (a7(A + m) + Zjv-i,iv-l) 
(A + m)2 — n{xp + Zjv^iv 


( 20 ) 


i-1 

where Ijj = 7 ^ X) Let us dehne the ith column of (20) as Cj. Let X = lu for any 

Zc=l 


non-zero 00 &TZ and apply the following elementary column operations: Ci — —t Ci 


and c ^ 


—cn —>■ Cj for 1 < j < iV. Thus, the matrix (20) is reduced to a matrix F 


aN — jry 


11 

















given by 


r 2 / \2 2 

—m — n{x) — u 

0 

0 

0 

+ m) 

/2,1 

/ 

0 

0 

a^~‘^ {a'y{iuj + m) + ^2,2) 

h,i 

h ,2 

/ 

0 

a^~'^ + m) + ^3^3) 


fN-1,2 

fN-1,3 

/ 

a {a'^poj + m) + In-i,n-i) 

/ni 

In ,2 

In ,3 

• • • fN,N-l 

{loj + m)^ — n(x)^ -|- In,n 


where / = — n{xY — xP' — jm + iuj{2m — 7 ) and fij (1 < i < iV, 1 < j < i) is a 

complex number. As a; 7 ^ 0 and —mP — n{x)‘^ — is a negative real number, it can be 
found that the columns of the matrix F are linearly independent, that is, for a vector 
[ui, ^ 2 , • ■ ■ , uat] e 7^^, the equality 


^^1 


F 


V2 


= 0 , 


( 22 ) 


vn 


is true if and only if r>i = r>2 = ■ ■ ■ = uat = 0. Thus, the matrix A„(A)Ai(A) — n{x)^lN 
has full rank, which leads to det (A„(A)A/(A) — u^xYIn) p 0 when a; 7 ^ 0. Consequently, 
det(zC(;/ — Aj^) 7 ^ 0 V non-zero 00 eTZ. Following Lemma we obtain the theorem. 

In the rest of paper, we shall numerically analyze stability and entanglement perfor¬ 
mance of the A^-NOPA coherent feedback system. To this end, we now set y = 1, that 
is, 7 = 7,.. Based on Theorem]^ with the help of Mathematica, we get Fig. [^that plots 
the values of the stability threshold Xth of our A7-NOPA systems {2 < N < 20) in the 
absence of losses (black circles), with transmission losses only (blue crosses) and with both 
transmission and amplification losses (red plus signs). The values of Xth of the A7-NOPA 
systems {2 < N < 6) are listed in Table It is indicated that the value of the stability 
threshold Xth decreases as more NOPAs are added to the network. The rate of decrease 
becomes smaller as the number of NOPAs grows. Moreover, existence of amplihcation 
and transmission losses broadens the range of x over which stability is guaranteed. Notice 
from the hgure that the effect of transmission and amplihcation losses on Xth diminishes 
for higher values of N. 

Note that when values of all system parameters except for x are given, we can also use 
the mussv function in the Robust Control Toolbox of MATLAB to estimate the stability 
threshold. Details are given in Appendix 1 . 


4 Entanglement 

In this section, entanglement performances are compared among the systems with number 
of NOPAs varying from 2 to 6 , in the absence of time delays. First, we study the EPR 
entanglement, namely, the two-mode squeezing, between the two outgoing helds when sys¬ 
tems are in an ideal case where no losses are present. After that, the effect of transmission 
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Figure 5: Values of Xth of iV-NOPA systems {2 < N < 20) in the absence of losses 
(a = 1, K = 0) (black circles), with transmission losses only (a = 10 , k = 0) (blue 

crosses) and with both transmission and amplihcation losses (a = 10 , k = j ^th) 

(red plus signs), with y = I and d = 1. 


Table 1: Values of stability threshold Xth of the V-NOPA coherent feedback system (2 < 
N < 6) without losses, with transmission losses only and with both transmission and 
amplihcation losses, y = I and d = 1. 


N 

Xth 

(a = 1, K = 0) 

Xth 

- 0.01 

(a = 10 ^-1 , K = 0) 

Xth 

2 

0.4142 

0.4209 

0.4363 

3 

0.2679 

0.2715 

0.2808 

4 

0.1989 

0.2013 

0.2080 

5 

0.1583 

0.1602 

0.1655 

6 

0.1316 

0.1331 

0.1375 
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and amplification losses on the two-mode sqneezing is taken into acconnt. In the end, 
we analyze the entanglement between pairs of optical cavity modes in the system nsing 
logarithmic negativity as an entanglement measnre. 


4.1 End-to-end EPR entanglement 


Strong correlation between qnadratnre-phase amplitndes of two helds is a manifestation 
of EPR entanglement [S]. The EPR entanglement between two continnons-mode helds 
is qnantihed in freqnency domain. To this end, we dehne the Fonrier transforms of 
^out,b{t) and iout,a{t) in Eq. as ^and = 

f^^out,a(i)e~'‘^^dt. The two-mode amplitnde spectrnm V+(zlj) and the two-mode 
phase spectrnm V-(zuj) are dehned as [8] 


((SL,a(*^) + = V+(zu;)6(cv - a;'), 

= R_(zu;)J(a; - u;')- 

Fnrthermore, based on (|^ and (|^, the two-mode sqneezing spectra are obtained by [2S1I2E] 


R+M = Tr[ifi(za;)*ili(za;)], (23) 

V.(zu)= TrlH2(zu)*H2(zu)], (24) 

where = [101 0]il(zu;), H 2 {iuj) = [010 —l]il(zci;), and H{ioj) = {lujl — ^v)~^ Bn+ 

Hat is the transfer fnnction of the iV-NOPA system. Note that V±{iu:) > 0 for all uj. 

Dehne the two-mode sqneezing spectrnm V{zoo) as 


R(zci;) = V^iiuj) + V_{iuj). (25) 

The helds ^out,a and ^outfi in the system shown in Fig. are EPR entangled at the freqnency 
oj rad/s if 39a, 9b € (—tt, tt] snch that V{toj, 9a, 9b) satishes the snm criterion [9], 

V{iu,9a,9b) = V+{iuj,9a,9b) + V.{iu,9a,9b) < A. (26) 


Perfect two-mode sqneezing has the featnre that V{zoo, 9a, 9b) = V±{zcj, 9a, 9b) = 0 [S]. 

Of conrse, perfect sqneezing cannot be achieved in practice. Therefore, one aims instead 
to have a small valne V±{zuj, 9a, 9b) over a wide freqnency range [9]. According to [I2l 1^ . 
V±{iuj, 9a, 9b) ~ V±(0, 9a, 9b) holds at low freqnencies. Thanks to this, we can simply focns 
on the two-mode spectra V (0, 9a, 9b) and Ri-(0, 9a, 9b) at a; = 0 for the rest of the paper. 

In this paper, plots of the two-mode spectra are presented in dB nnit, that is, Ri-(za;)(dB) = 
10 log;^o V±{zuj) and V (za;)(dB) = 10 log^g ^ f^is case, perfect EPR entanglement at 
freqnency uj corresponds to V±{zuj) = —oo (dB). Better two-mode sqneezing is indicated 
by a more negative valne of R(zC(;)(dB). 


4.1.1 An ideal case. 

Now we examine the two-mode spectra when the A-NOPA system is lossless. First, we 
aim to hnd valnes of 9a and 9b at which the cost fnnction V (0, 9a, 9b) at a; = 0 is minimized. 
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Based on (23), (24) and via Mathematica, we obtain Table ?? which presents the formulas 
of the two-mode squeezing spectra of iV-NOPA systems in the ideal case. 


For any value of x in the interval (0, Xth), we have — 1 -|-x^ < 0 for the 2-NOPA system, 
3 — lOa:^ -|- 30:“^ > 0 for the 3-NOPA system, — 1 -|- 7x‘^ — Tx"^ -|- x® < 0 for the 4-NOPA 
system, 5 — 60x^ -|- 126x'^ — 60x® + 5x® > 0 for the 5-NOPA system, and —3 -|- 55x^ — 198x^ -|- 
198x® — 55x® -|- 3x^° < 0 for the 6-NOPA system. Recall that the values of the stability 
threshold Xth are listed in Table [Tj Thus, for systems with an even number of NOPAs, the 
best two-mode squeezing is obtained if 9a + Ob = 0 or 9a = 9b = 7i; for systems with an odd 
number of NOPAs, the best two-mode squeezing is achieved when -|- 6'b| = vr. In this 
regard, we set 


J (0, 0), for N is even, 
(tt, 0), for N is odd. 


(27) 


In the rest of paper, V{iuj) and V±{iu) are dehned as the two-mode spectra at the hxed 
values of 9a and 9b as given in ( [^ . 

Now we check the two-mode spectra V±(0, k) as a function of /c at a; = 0, where x = kxth 
as the value of k varies from 0.5 to 1. As shown in Fig. the two-mode squeezing spectra 
decrease as the value of x approaches the stability threshold. Moreover, the rates of the 
decreases are similar. 

One of our interests is the power consumption of the systems to generate the same level 
of EPR entanglement, say, R(0) = —25 dB. We denote the corresponding x as x _25 dB- 
Here the power of pump beam employed by each NOPA is Hence, the total pump 

power of an A^-NOPA system is Nx^'y^. Since is a fixed reference, to compare pump 
consumption between different values of N, it is enough to consider only the quantity Nx"^. 
As the third column of Table indicates, a system with more NOPAs consumes less power 
to yield the same degree of two-mode squeezing. 


4.1.2 Effects of losses. 

The effect of losses on the two-mode squeezing of the systems is indicated in Table All 
the systems have the same value of V (0) {V (0) = —25 dB) when losses are neglected. EPR 


Table 2: Two mode squeezing spectra of the A^-NOPA coherent feedback system {2 < N < 


6) without losses, y = 1. 


N 

n(o) 

2 

r) (l+x^)'^ + 16x^( — l+x"‘y+8x{l+x"‘y( — l+x^) COs{9a+dh) 

(1 — 

3 

r, (l+a;^)'^+4x^(3-10x'-'-|-3a:"‘)^+4a;(l+x^)^(3-10a;^+3x'*) cos(6»a+6»6) 

(—H-15x^ —15x'^+x®)^ 

4 

r) (l+x'‘‘')^+64x'‘‘(—l+7x'‘‘—7x'^+x'^y‘‘+16x{l+x'^y{—l+7x‘'‘—7x‘^+x*^) cos{da+9b) 

^ (l-28a:2+70x4-28a:6+a;«)2 

5 

D(l-|-a;^)^'J+43:^(5-603;^+126a;'‘-60x‘>+5x«)^+4a:(l+a;^)^(5-603;^ + 126a;"‘-60a:‘’+5a;»)cos(6»a+6»i,) 

^ (-l+45x2-210a:'‘+210a;6-45x8+3;lO)2 

6 

o(l+a:^)^^+16x2(-3+55a;2-1983:'‘+198a;'>-55x»+3a:^^)2+8x(l-|-a;2)0(_3+55a;^_198x'‘+198a;t’-55a;«+3a:^^)cos(6»a+6»6) 

(1—66x^+495x'^—924x®+495x^— 
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Figure 6: Log-log plots of 1/(0, fc) (dB) with respect to iV-NOPA systems {2 < N < 6 ) 
with k varying from 0.5 to 1, x = kxth, y = I, a = 1 and k = 0. 


entanglement of each system is degraded by around 15 dB under the effect of transmission 
losses, and the reduction is more than 20 dB if both transmission and amplihcation losses 
are present. Generally, a system employing more NOPAs provides a slight improvement 
in EPR entanglement when transmission losses are present, in the absence of amplihcation 
losses. This merit disappears in the presence of both transmission and amplihcation losses. 
In this case, the system with more NOPAs yields less EPR entanglement. 


Table 3: Power consumptions {Nx'^), values of V±(0), and values 1/(0) of A^-NOPA systems 
(2 < N < 6 ) under ehect of losses, with x = X -25 dB, 2/ = 1 and d = 1. 


N 

X-2b dB 
(a = 1, 

K = 0) 

^^-25 dB 

F±(0) 

-0.01 

[a = 10 N-i , 

K = 0) 

F(0) 

-0.01 

(a = 10 ^-1 , 

K = 0) 

F±(0) 

-0.01 

[a = 10 ^-1 , 

_ 7 3xl06 \ 

F(0) 

-0.01 

(a = 10 ^-1 , 

_ 7 3xl06 \ N 

[o.6xV2)^> 

[o.6xV2)^> 

2 

0.3978 

0.3165 

-13.3150 

-10.3047 

-7.5838 

-4.5735 

3 

0.2579 

0.1995 

-13.3286 

-10.3183 

-7.4114 

-4.4011 

4 

0.1916 

0.1468 

-13.3302 

-10.3199 

-7.3510 

-4.3407 

5 

0.1526 

0.1164 

-13.3295 

-10.3192 

-7.3236 

-4.3133 

6 

0.1269 

0.0966 

-13.3306 

-10.3203 

-7.3078 

-4.2975 


Now we compare the two-mode squeezing levels when the systems are consuming the 
same total pump power. From now on, we use xn to denote x of the system with N 
NOPAs. In this case, we set x% = 0.13, hence Xi = {^JQ/^)x% {i = 2, 3,4, 5). The two-mode 
squeezing spectra in the (ideal) lossless case, in the presence of transmission losses only. 
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Table 4: lApproximate optimal two-mode squeezing under the effect of transmission losses 

and the corresponding power consumption for the A^-NOPA systems {2 < N < 6) with 
- 0.01 

y = 1, d = 1, a = 10 ^-1 and k = 0 . _ 


N 

^opt 

TV 

‘^opt 

14 (0) 

no) 

2 

0.4074 

0.3319 

-13.3683 

-10.3580 

3 

0.2644 

0.2097 

-13.3928 

-10.3825 

4 

0.1965 

0.1544 

-13.3991 

-10.3888 

5 

0.1565 

0.1225 

-13.4018 

-10.3915 

6 

0.1302 

0.1017 

-13.4033 

-10.3930 


Table 5: Approximate optimal two-mode squeezing under effects of both transmission and 
amplification losses and the corresponding power consumption for the A^-NOPA systems 
{2 < N < 6) with y = 1, d = 1, a = 10“*^ and n = x- 


N 

^opt 

N 

‘^opt 

>4(0) 

V(0) 

2 

0.3770 

0.2843 

-7.6545 

-4.6442 

3 

0.2435 

0.1779 

-7.4982 

-4.4879 

4 

0.1805 

0.1303 

-7.4435 

-4.4332 

5 

0.1438 

0.1034 

-7.4182 

-4.4079 

6 

0.1195 

0.0857 

-7.4044 

-4.3941 


as well as under effect of both transmission and amplification losses are plotted in Fig. 
With the same total power, a system consisting of more NOPAs yields stronger EPR 
entanglement, except that when both transmission and amplification losses are present, 
the 6-NOPA system has a slightly lower degree of EPR entanglement than the 5-NOPA 
one. Moreover, the EPR entanglement of the system with less NOPAs has smaller change 
under the influence of losses. 

To find the approximate value of Xopt at which the system achieves the highest degree 
of two-mode squeezing at ca = 0 under the effect of losses, we pick the smallest one among 
the values of V±(0) corresponding to a thousand samples of x evenly spread through the 
range [0.001, l]xth- Tables and illustrate the values of Xopt, the corresponding two¬ 
mode squeezing degrees, and the total power consumptions of the A^-NOPA systems in the 
scenarios with only transmission losses as well as when both transmission and amplification 
losses are present, respectively. As the tables indicate, the best two-mode squeezing degrees 
of all the systems are similar, however, the system with more NOPAs consumes less total 
pump power. For instance, a 6-NOPA system needs less than a third of power used by the 
dual-NOPA system. Thus, the system should employ more NOPAs in the presence of losses 
for efficient use of pump power, while only losing a small amount of EPR entanglement. 
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Figure 7: Log-log plots of V±{ioj) (left) and V+iiuj) + V-{iuj) (right) of iV-NOPA systems 

(2 < N < 6) without losses (top, a = 1, k = 0), with transmission losses only (middle, a = 
- 0.01 - 0.01 
10 ^-1 , K = 0) and with both transmission and amplihcation losses (bottom, a = 10 , 

K = x), under the same total pump power, with xq = 0.13, Xi = {^/ 6 Ji)xQ 

{i = {2, 3,4, 5}), y = 1 and d = 1. 
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4.2 Entanglement of two-mode Gaussian states 


In this sub-section, we study the entanglement of two-mode Gaussian states with respect 
to the cavity mode operators when the iV-NOPA system is lossless. To this end, we hrst 
calculate the covariance matrix P^it) and steady-state covariance matrix P/v of fhe 2N- 
mode Gaussian state of the system. Based on Eq. (|^, 0 and ( [l0| ), we employ the Matlab 
functions ode^S with the sampling time 10“^° sec and lyap to compute PAr(t) and P^r, 
respectively. Here, we take Pq = I corresponding to that the system starts in a vacuum 
state. 


The covariance matrix PN,ai,b{t) of the modes ai and bj (i,j = {1,2, ■■■ ,N}) is a 
corresponding 4x4 sub-matrix of P/v(^)- For instance, the covariance matrix of Oi and 
in a 6-NOPA system is 


(f 6 )Ll 

(n)l,2 

(-^6)1,23 

(P6)l,24 

{Pnhi 

(-^5)2,2 

{ P &) 2,23 

(-^3)2,24 

(-^6)23,1 

(-^6)23,2 

{ P &) 23,23 

(-^6)23,24 

(-^6)24,1 

(-^6)24,2 

(-^6)24,23 

(^6)24,24 


(28) 


Entanglement of two-mode Gaussian states with corresponding covariance matrix PN,ai,bj {t) 
is measured by the logarithmic negativity Ej^^ai,b (t) [TOl |23]. Write PN,ai,b {t) in a 2 x 2 
block matrix form given by 


PN,ai,bj{t) 


PN,ai,bj,lif) PN,ai,bj,2if) 
PN,a„bj,2it)'^ PN,a„bj,3it) \ ’ 


(29) 


where ^k{t) {k = (1, 2, 3}) is a 2 X 2 matrix. Dehne 


AN,ai,bAt) = det{PN,ai,bj,l{t)) + det(Pjv^ai^,3(t)) - 2 det(Pjv^a. 2 (t)), 


J^N,ai,bj(t) - 


^N,ai,bj{t) - J AN,ai,bj{ty - A det{PN,ai,bj (t)) 


Then E]s[^ai,bj{t) is a nonnegative real number given by 

EN,ai,bj{t) = max[0, - log 2 l'N,ai,bj{t)\- 


(30) 

(31) 

(32) 


EN,ai,bj{t) = 0 represents that the modes and bj are separable at time t, that is, there 
is no entanglement between the modes. Strong entanglement between modes a* and bj is 
represented by a high value of E^^ai,bj{t). 


What is of interest to us are the time evolution and steady-state values of logarithmic 
negativities PAr,a,,b,(t), pAr,ai,fei+i(t), EN,ai+ubi(t), P7v,ai,bjv(^) as well as EN,aM,bi(t) of an N- 

NOPA {2 < N < 6) coherent feedback network. Besides, the logarithmic negativity of Oc = 
N N n 

tti and be = Y also looked into. The reason is that, as (2) indicates, when 
2=1 2=1 

losses and delays are neglected, the outputs ^out,a and ^out,b contain Oc and be, respectively. 
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Table 6: The steady-state logarithmic negativities of iV-NOPA systems {2 < N < 6) 
in the absence of losses and delays, under the same total pump power, with xq = 0.13, 
Xi = {\/Q/i)x6 {i = {2, 3,4, 5}), y = 1 a = 1, and k = 0. 


N 

E/2,ac,bc 

A'2,ai,f)2 

A'2,a2,fei 

-^ 2 ,ai ,61 

-£"2,02,62 



2 

0 

0 

0.4850 

0.1921 

0.1921 



N 


.^3,ai,f)2 

.^^3,02,63 

-£' 3 , 03,61 




3 

0.0561 

0 

0 

0.2865 





A'3,a2,&i 

A'3,a3,b2 

^"3, 01,63 






0.3645 

0.3645 

0 







A'3,a2,b2 

-^3,03,63 






0.1144 

0.1144 

0.1144 





N 

E^^acbc 

-^4,ai ,&2 

^<4,02,63 

-^4,a3,64 

-^4,a4,&i 



4 

0 

0 

0 

0 

0.1843 




-^4,a2 

774,a3,f)2 

-^4,a4,fe3 

-^4,ai ,64 





0.2803 

0.3021 

0.2803 

0 





-^4,ai 

A'4,a2,f)2 

-^4,03,63 

-^4,a4,64 





0.0722 

0.0722 

0.0722 

0.0722 




N 

E5^ac,bc 

A'5,ai,f)2 

-^5,02,63 

-£"5,03,64 

-£" 5 , 04,65 

-£"5,05,61 


5 

0.0223 

0 

0 

0 

0 

0.1207 



E^^a2,b\ 

.^5,a3,f)2 

^'5,04,63 

-£"5,05,64 

-£" 5 , 01,65 




0.2134 

0.2500 

0.2500 

0.2134 

0 





A'5,a2,f)2 

-^5,03,63 

-^5,a4,b4 

-£' 5 , 05,65 




0.0451 

0.0451 

0.0451 

0.0451 

0.0451 



N 

Ee^acbc 

-^6,ai,b2 

-^6,02,63 

-£" 6 , 03,64 

-£" 6 , 04,65 

-£"6,05,66 

-£"6,06,61 

6 

0 

0 

0 

0 

0 

0 

0.0767 


A'6,a2,bi 

A'6,a3,f)2 

^'6,04,63 

-£" 6 , 05,64 

-£" 6 , 06,65 

-£"6,01,66 



0.1552 

0.2033 

0.2195 

0.2033 

0.1552 

0 



-^ 6 ,ai , 6 i 

A'6,a2,f)2 

-£'6,03,63 

-^ 6 ,a 4,64 

-£" 6 , 05,65 

-£"6,06,66 



0.0260 

0.0260 

0.0260 

0.0260 

0.0260 

0.0260 
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Figure 8: Time evolution of E^^acbc (2 < iV < 6), in the absence of losses and delays, 
under the same total pump power, with xq = 0.13, Xi = {^JQ/i)xQ {i = {2, 3,4,5}), y = 1, 
a = 1, and k = 0. 


Notice that Oc and be can be viewed as collective single mode annihilation operators as 
they satisfy the commutation relations [oc, a*] = 1, [^o^c] = 1) ^c] = 0 and [oc, &*] = 0. 

The covariance matrix of Oc and he is PN,a^,be(t) = ^nPn with = ^[J ■ ■ ■ /]. 

In the absence of losses and delays. Table j^shows the values of steady-state logarithmic 
negativities of the iV-NOPA systems. Fig. [^indicates the time evolution of i?jv,ac,6c) 

Fig. 1^ plots the evolution of logarithmic negativities of a 6-NOPA coherent feedback net¬ 
work. As indicated, and 6j+i remain separable for all time, and the same happens to oi 
and h]^. At steady state, entanglement exists between modes bi and Oj+i, and bi, as well 
as Oat and bi. In particular, it can be observed that internal entanglement synchronization 
occurs at steady state, that is, the degree of entanglement between the oscillator modes 
tti and bi in the cavity of each NOPA in the A^-NOPA coherent feedback network is the 
same. For systems with an odd number of NOPAs, there is slight entanglement between 
the collective modes Oc and be, while in systems containing an even number of NOPAs, Oc 
and be are entangled at the beginning for a very short time. After that, the entanglement 
rapidly vanishes. Moreover, with the same total pump power, entanglement of two-mode 
Gaussian states in the system with more NOPAs is weaker. It is an interesting result that 
even though the two-mode entanglement of the internal cavity modes does not improve 
for systems carrying more NOPAs, its EPR entanglement between the two outgoing helds 
does improve. 
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Figure 9: Time evolution of logarithmic negativities of the 6-NOPA system in the absence 
of losses and delays, under the same total pump power, with xg = 0.13, y = 1, a = 1, and 

K = 0. 
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5 Effect of Time Delays 


In this section, we investigate stability and entanglement of the iV-NOPA systems in the 
presence of losses and time delays. For a d km transmission distance, the time delay r of 
each path between two neighbouring NOPAs is r = 3 xio- 5 (jv-i) • check stability of our 
time-delayed A^-NOPA systems, we employ the DDE-BIFTOOL toolbox [271 EH], which 
is a Matlab package used to plot the eigenvalues of a linear delay differential system. A 
system is stable if all real parts of the eigenvalues are negative. Based on the above fact, 
we hnd that stability of the A^-NOPA systems for N up to six is guaranteed in the case 
where the systems are given the same total pump power, xq = 0.13 and both losses and 
delays are present. 


As a linear quantum system, the A^-NOPA network with time delays can be built 
in Matlab via commands connect and delayss in the Matlab Control System Toolbox. 
The non-rational transfer functions (due to the time delays) Hi{s) and 7^2(s) in (23) and 
(24) can be numerically computed with the built-in Matlab frequency response command 
freqresp. Therefore, the two mode squeezing spectra V±{iuj) are obtained via (23) and (24). 
The effect of time delays on EPR entanglement between the outgoing helds of our A^-NOPA 
system is indicated in Fig. IT and Fig. H, where all the systems are given the same total 
pump power and undergoing both transmission and amplihcation losses. Compared with 
the two-mode squeezing of the systems in the absence of delays, the presence of time delays 
reduces the bandwidth over which the EPR entanglement exists, but does not impact the 
EPR entanglement degrees at low frequencies. The phenomenon of the sharp peaks and 
dips at high frequencies is a common feature of the frequency response of systems under 
the effect of internal time delays, see, e.g, [29] p. 182]. In our case, the bandwidth of EPR 
entanglement under influence of time delays is similar for all the systems with a different 
N. 


6 Conclusion 

This paper has studied the stability condition and entanglement performance of an N- 
NOPA coherent feedback network with N up to six, where the NOPAs are evenly dis¬ 
tributed in a line between two distant parties, Alice and Bob. The system undergoes 
transmission losses, amplihcation losses and time delays. Moreover, two adjustable phase 
shifts 6 a and 6 h are placed at Alice and Bob for achieving the best two-mode squeezing 
between the two outgoing helds by selecting appropriate quadratures of the output helds. 

In the absence of time delays, we have derived a necessary and sufficient stability 
condition with the aid of /x-analysis method from control theory [22| by regarding x, 
the parameter related to the amplitude of pump beam, as an uncertainty. We have shown 
that, the value of stability threshold Xth is the smallest positive root of the polynomial 
det (Ajv(a;)). It is observed that the existence of losses broadens the range of x over which 
stability is guaranteed, and the value of Xth decreases as more NOPAs are added to the 
system. 
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Ns2, with time delays. 




Figure 10: Log-log plots of V±{iu) (left) and V{iu) (right) with respect to a 2-NOPA 
system (top), a 3-NOPA system (middle) and a 4-NOPA system (bottom), without time 
delays (blue solid line) and with time delays (magenta dashed line), under the same total 
pump power, in the presence of losses, with x% = 0.13, Xi = {^JQ/i)xQ {i = {2,3,4}), 

« = 10^, K = x,T= 3^10-Wd ’ 1/ = 1 and d = 1. 
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Figure 11: Log-log plots of V±{iuj) (left) and V{iu) (right) with respect to a 5-NOPA 
system (top) and a 6-NOPA system (bottom), without time delays (blue solid line) and 
with time delays (magenta dashed line), under the same total pump power, in the presence 

of losses, with xg = 0.13, Xg = a = 10^, k = x, r = 3 xio- 5 (jv_i) , 

y = 1 and d = 1. 
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Strong EPR entanglement is represented by strong attenuation of two-mode output 
squeezing spectra below the sum criterion. In the ideal case, we have found the values 
of 6a and 6h at which the system achieves the best two-mode squeezing. Moreover, the 
two-mode squeezing increases rapidly as the value of x approaches the stability threshold 

Xth- 

We have compared the two-mode squeezing generated by systems with different numbers 
of NOPAs. It is shown that, to achieve the same squeezing level in the ideal case, the system 
employing more NOPAs requires less total pump power. When losses are present, all the 
systems have a large and similar decrease in EPR entanglement. Given the same total pump 
power, the system carrying more NOPAs has improvement in the two-mode squeezing in 
the ideal case and when only transmission losses are present. However, this is no longer 
assured when amplihcation losses are also taken into account. Furthermore, the best two¬ 
mode squeezing degrees of the systems with losses are similar. However, the system with 
more NOPAs requires less pump power to achieve the best two-mode squeezing. 

We have also investigated the entanglement of two-mode Gaussian states of the internal 
cavity modes. Steady-state values and time evolution of logarithmic negativities have been 
studied. It is shown that entanglement exists between the modes a* and 6*, bi and Oj+i as 
well as bi and a^- Moreover, we have observed an internal entanglement synchronization 
that occurs between the modes ai and bi for i = 2, 3,4,5, 6 at steady state. In the ideal 
case, given the same pump power, though the system with more NOPAs has improved 
two-mode squeezing in the output helds, it does not have better internal entanglement 
between cavity modes as measured by logarithmic negativity. 

Stability and entanglement under the effect of time delays has been studied as well. 
With time delays, stability is checked with the DDE-BIFTOOL toolbox [271 l28] . It is 
shown that, with transmission and amplihcation losses, time delays narrow the bandwidth 
over which the EPR entanglement exists. 

This work gives several qualitative Endings on entanglement in the A^-NOPA network 
when N > 2, which have not been quantitatively analyzed and are topics suitable for 
future investigations. It is observed that the two-mode squeezing spectra decreases as the 
value of X approaches Xth when the system is lossless, and there exists an optimal value of 
X for the two-mode squeezing when the system is under losses. These phenomena are not 
surprising and to some extent predictable as they were proved for the 2-NOPA system in 
our previous work, see Theorem 2 and Theorem 3 in [T2]. In addition, future work is still 
required to quantitatively analyze the observation that entanglement between collective 
modes behaves differently between systems with an even and the odd number of NOPAs; 
also, adding more NOPAs into the system improves the end-to-end entanglement between 
continuous-mode output helds but not the entanglement of internal two-mode Gaussian 
states when the system is lossless and consumes the same pump power. Moreover, it would 
be of interest to quantitatively analyze the behaviour of the linear coherent feedback chain 
in the asymptotic limit of A^ —)■ oo. 
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Appendix 1 


When values of all parameters except for x of an iV-NOPA coherent feedback system are 
given, we can employ the mussv function in the Robust Control Toolbox of MATLAB 
to approximate the stability threshold. The approximated stability threshold, say Xthi 
approaches the stability threshold Xth from the left. That is, the system is robustly stable 
when the value of x belongs to the range {Q,Xth\ and Xth approximates the threshold value 
Xth- The value of Xth is found by the following bisection algorithm. 

Step 1. Start from Xth = 1. If the system is not stable over the range x G (0,a;t/i], set 
Xh = Xth and xi = 0; otherwise stop. 

Step 2. Set Xth = If the system is not stable over the range x G set 

Xh = Xth') otherwise set xi = Xth- 

Step 3. If the value oi Xh — xi > e for a prespecihed error tolerance e > 0 (here, we take 
e = 10“^°), go back to Step 2; otherwise set Xth = xi, stop. 
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